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GENERATION OF PSEUDO-RANDOM NUMBERS 

I. INTRODUCTION 


The advent of high-speed digital computers has made it possible to use simula- 
tion techniques incorporating probabilistic features. These simulations are generally 
referred to as Monte Carlo simulations and are resorted to whenever the systems being 
studied are not amenable to deterministic analytical methods or where direct experi- 
mentation is not feasible. An integral part of these simulations is the use of random 
numbers having a certain specified distribution characteristic of the process being 
studied. Since these simulations require a vast amount of calculations, speed is a 
vitally important factor. It was soon recognized that it was intractable to fill the 
computer's memory with large aggregates of random numbers. Because of this, 
computer algorithms were developed which allowed random numbers to be generated 
on-line. Since the generation of random numbers by such numerical algorithms is 
somewhat a contradiction in terms, they are often called "pseudo -random" numbers. 
There are essentially two problems encountered in generating such random numbers. 
One is that the generated random numbers are not representative of the desired 
distribution, and the other is that they are not statistically random, i.e., that there 
exist correlations in the generated numbers. The latter problem is the more serious 
one as evidenced by the considerable attention it has been given in the literature. 

This report presents several random number generators which have been found 
particularly useful in aerospace engineering applications. APL computer programs are 
also listed in Appendix B for most of the generators. 


II. UNIFORM RANDOM NUMBER GENERATORS 


The basic element of all Monte Carlo simulations is the uniform random number 
generator. Once uniform random numbers are available, all other desired distributions 
can be obtained either by use of the probability integral transformation or by applying 
some known relationship between the desired distribution to be generated and the 
uniform distribution.! One of the early methods used to generate uniform random 
numbers was the mid-square method originally proposed by John V. Neumann. In 
practice, one selects an arbitrary K-digit number, squares it, and then selects the K 
middle digits as the new random number. The process is repeated using this new 
random number. The drawback of this method is that it can produce a zero random 
number at unpredictable times and, thus, the process terminates. Consequently, 
this method was abandoned quite early in favor of the so-called congruential method 
first proposed by D. H. Lehmer in 1949 [1]. Accordingly, the random number 
generator takes the form 

17 Often it is convenient to generate a random number from a specified distribution 
by employing its relationship to the normal distribution. However, one must first 
generate the normal random number (s) as a function of uniform random numbers 
and then proceed to the desired distribution. 


( 1 ) 


E (a X. + b) Mod M , 


where the multiplier a, increment b, and modulus M are integers.. The starting value 
X^ is called the "seed" of the random number generator. Since the congruential 

relationship (1) is cyclical, the sequence of random numbers will repeat after a cer- 
tain period. Generators in which b = 0 are called multiplicative; otherwise, they are 
called mixed. The statistical behavior of the generated random numbers is predomi- 
nantly governed by the choice of the multiplier a and the modulus M. Therefore, the 
most widely used generators are of the multiplicative type (b = 0) . Many empirical 
and theoretical tests have been developed to assess the "goodness" of a random 
number generator. One of the most popiolar of these tests today is the lattice test 
which determines the lattice structure of the random number generator by comparing 
successive n-tuples (x. , ..., x.^^_^) and (x.^^, x.^2’ •••’ According 

to this test an acceptable generator can be obtained by selecting the constants a, b, 
and M to achieve a nearly hyper-cubic lattice structure, i.e. , the ratio of cell sides 
should be close to unity [2]. In practice, n is usually less than or equal to five [3]. 
Especially useful are congruential generators for which the modulus M is a prime 
number and hence, are called prime modiolus generators. If the multiplier a is selected 
to be a primitive root modulo M , the generated random number sequence attains its 
maximum period P = M-1 and all possible value from 0 to M-1 will be generated. 

Two useful uniform random number generators of this type which have very 
satisfactory lattice structures are: 


a = 7^ = 16,807 , b = 0 , M = 2^^-l ^ 2,147,483,647 , (2) 

with 5- dimensional cell side ratios 1:7.60:3.39:2.09:1.67, and 


a = 7^02,479 jv,o(j(2^^-l) = 29,903,947 , b - 0 , and 

M = 2^^-l = 2,147,483,647 , (3) 


with 5-dimensional cell side ratios: 1:1.04:1.30:1.22:1.09 [4], 


III. RANDOM NUMBERS FROM CONTINUOUS DISTRIBUTIONS 


Many of the generators for continuous distributions are obtained by a direct 
application of the probability integral transformation [5]. For a given uniform random 
number u between zero and one, a random number x having the desired distribution 
F(x) is obtained by solving the equation u = F(x) for x. Since this process requires 

the determination of the inverse cumulative distribution function F ^(x) , its prac- 
ticality depends upon the availability of explicit expressions or convenient approxima- 
tions for this inverse cumulative distribution function. We now discuss methods on 
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how to generate random numbers from continuous distributions which appear frequently 
in aerospace engineering simulations. 


A. The Normal Distribution 

The most common distribution is the normal or Gaussian distribution with den- 
sity function given by 


f(x;y ,a^) 



-“ < X, y < <» and a > 0 . (4) 


Consequently, a variety of methods have been devised which can be used to generate 
normal random numbers. Of these, the following three were found to be satisfactory 
and practical in terms of accuracy and computer run time. 


1. Hastings Ap proximation 

This method invokes the probability integral transformation in a slight variation 
in that it employs the complement of the cumulative distribution Q(x) = l-F(x). The 
reason for this is that smtable rational approximations for Q(x) were derived by 
Hastings [6]. The most accurate of Hastings approximations is given as: 


X = t 


C + C, t + C 
o 1 

1 + dj^t + d2t^ 



+ dst 


3 


+ e(p) 


(5) 


where 


C = 2.515517 
o 

= 0. 802853 
= 0.010328 

d^ = 1.432788, t = /21np 
dg = 0.189269, 0 < p <0.5 

= 0. 001308 , and | e(p) | < 4. 5 x 10 


X is the desired Normal random number and p is a uniform random number. 

2. Box-Muelle r or ^Tolar^^ Method 

This method generates a pair of normal random variables using a pair of uni- 
form random numbers as follows : 
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Let and U 2 be independent uniform random variables and define 


1/2 

= (-21n u^) cos 27 t , 


1/2 

%2 = (-21n Uj^) sin 2 it U 2 


( 6 ) 


Then and X 2 are two independent normal random variables with zero mean and unit 
variance. To see this, we establish the inverse relationships 

^ 2 ^ 2, 

(Xf + X2 ) 

2 

u^ = e 


and 


(7) 


u„ = ^ tan ^ f — 
2tt I Xj^ 


It follows then that the joint probability density function of Xj^,X 2 is 


/ 2 ^ 2, 
(Xf + X2 ) 


f (x^,X2) 2 tt ® 


= f(x^) f(X 2 ) 


( 8 ) 


and, thus, the desired conclusions, including the independence of x. and x„, are 
obtained [ 7] . 

3. Central Li mit Me'yiqd 

Let X., x„, ... , X be a sequence of n uniform random variables. Then 
j. z n 


y 


n 




(9) 


will be distributed asymptotically as a normal random variable with zero mean and unit 
variance. For n = 12, we see that (9) reduces to 

12 

5^12 = ® • 

i=l 
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The exact distribution of the standardized sum of n independent uniform random 
numbers can be easily derived using moment generating functions. Since the sum x 
of n independent uniform variables has moment generating function 


~ (1 - e'b” 

X 


( 10 ) 


its density function is given as 


f(x) = 


(n-1)! 


k=0 


[(x-k) u(x-k)] 


n-1 


( 11 ) 


where 


u(s) = 


0 

1 


s < 0 
s > 0 


and 0 X n 


For the standardized random variable y^ in equation (9), we find that 



where - y§rT < y < i/§n 

Density functions obtained by this method for n = 2, 4, 12, 20 are compared 
with the normal density fimction (dotted) in Figure 1. Also, a comparison of the 
cumulative distribution function for n = 12 with the cumulative normal distribution is 
given in Table 1 to four decimal places. 

The agreement between the two distributions is very good except in the tail 
areas. But a comparison of random numbers generated by the three methods revealed 
no significant statistical differences even for the tail areas. 

A comparison of computer CPU run times for each of these three methods to 
generate 1000 normal random numbers is as follows : 
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Method 

Time 

Hastings Approximation 

Box -Mueller 

Central Limit Method 

1.13 sec 
1.02 sec 
1.08 sec 





Figure 1. Normal density function (dotted) and standardized 
sum of uniforms (n = 2, 4, 12, 20). 

TABLE 1. CUMULATIVE NORMAL DISTRIBUTION 


X 

F(x) 


-4.0 

3.167 ^ 10'^ 

8,526 X 10'® 

-3. 5 

2.326 ^ 10'^ 

1.212 X 10“'^ 

-3.0 

1.350 N 10"^ 

1.007 X 10'^ 

-2. 5 

6.210 ^ 10'^ 

5. 579 X 10'^ 

-2.0 

2.275 X 10'^ 

2. 227 X 10'^ 

-1.5 

6.681 X 10"^ 

6.745 X 10'^ 

-1.0 

1.587 X 10'^ 

1.608 X 10'^ 

-0.5 

3.085 X 10"^ 

3,106 X 10'^ 

0.0 

5.000 X 10'^ 

5.000 X 10'^ 

0.5 

6.915 X 10'^ 

6. 894 X 10“^ 

1.0 

8.413 X 10'^ 

8.393 X 10*^ 

1.5 

9.332 V 10'^ 

9.326 X 10'^ 

2,0 

9.773 X 10'^ 

9. 777 X 10'^ 

2. 5 

9. 938 V 10'^ 

9. 944 X 10'^ 

3.0 

9.987 X 10'^ 

9. 990 X lo"^ 

3.5 

9. 998 X 10'^ 

9. 999 X 10’ ^ 

4.0 

1.000 

1.000 


♦Comparison of the cumulative normal distribution 
F(x) with the cumulative distribution function for 
the central limit method using equation (13) as the 
density function with n = 32. 






B . Log-Noririal Distribution 


It is often claimed that the log-normal distribution is as fundamental as the 
normal distribution and may be thought of as arising from the combination of random 
terms by a multiplicative process. The log-normal distribution has been applied in a 
wide variety of fields including social sciences , physical sciences , and engineering and 
its density function is given by 


f(x;y,a) = (2'na^-x.^) exp -1/2 


Log-normal random numbers x may be generated by the relationship x = e^, where y 
is a normal random number obtained by methods discussed in Section III, paragraph A. 




0 < X , a < • 

- CO C 1 1 <* CO 


(13) 


C . Weibull Distribution 

The Weibull distribution is perhaps the most popular distribution at the present 
time when dealing with problems of reliability and material fatigue. Its appeal stems 
from its mathematical tractability and for this reason is often preferred to the gamma 
distribution. The Weibull density function is given by 


f(x;a, g) = ag x^ ^ exp(-a x^) , 0<x<<» a, g>0 . (14) 

— 0 

The cumulative distribution function F(x) — 1 -e leads immediately to the inverse 
relationship 


= ( inu ) 


as the desired Weibull random generator. 


(15) 


D. Gamma Distribution 

The gamma distribution is another two parameter disti?7.bution which is also quite 
flexible in fitting a variety of random processes. It finds use in reliability analysis, 
meteorology, and aerospace engineering. Its density function is 


f(x) 



-OCX g- 1 

e X 


( 16 ) 


where 0 < x < <» and a,g > 0. 



lillllll llllllllll■llllllllllllll llllll■lllllllllllllllllllll 


For integer values of 3, the gamma distribution is often referred to as the Erlangian 
distribution after the Danish mathematician, A. K. Erlang, who introduced it in the 
theory of queues and Markov processes in 1917. Random numbers following the 
Erlangian distribution are generated by the formula 


X 



(19) 


E. Exponential Distribution 

The exponential distribution appears often in engineering applications because of 
its importance in reliability theory and queueing theory. The probability density 
function of the exponential distribution is 


f (x ; a) = a e 


•C(X 


( 20 ) 


where x,a > 0. In this case, an explicit expression for the inverse cumulative dis- 

- ax 


tribution function is easily obtained by solving the equation u — 1 
producing 


for X, 


X 


In u 

a 


( 21 ) 


as the desired exponential random number generator. 


F. Chi-Square Distribution 

2 

As a special case of the gamma distribution, the x -distribution is often used 
as a measure of goodness of fit of a specified distribution to observed frequencies. 

2 

For discrete variates , x provides a sensitive test of departure from the Poisson 
distribution. The Chi-Square density function with n degrees of freedom is given by 


f (x;n) 


X n-1 ^-x/2 
2 " r (n) 


^ > 0 and n = 1,2,3,... 


( 22 ) 


Because of its relationship to the normal distribution, Chi-Square random numbers 
may be generated by taking the sum of n squared Normal random numbers. 
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G . F Distribution 

The F distribution appears extensively in statistical hypothesis testing under 
normality theory and has the density function 


f(x;m;,n) 


riJ2 mJ2 
n m 



1 


r(n/2) r(m/2) (m + nx) 


(23) 


where x > 0 and m,n = 1,2,..., are degrees of freedom. In practice, it is generally 
found to be most convenient to make use of the F distribution's relationship to the 

2 

Chi-Square distribution in order to generate F random numbers. That is, if x (™) 

2 

and X (ii) two independent Chi-Square random variables with m and n degrees of 
freedom, respectively, then 


2 

Y = X (m)ym 
X^(n)/n 


(24) 


follows the F -distribution with m, n degrees of freedom. Thus, 


m 



i=m+l 


(25) 


is the desired F random number generator with m and n degrees of freedom and the 
X^, i = l,2,...,m+n are Normal random numbers. 


H. Beta Distribution 

To generate random numbers from the Beta distribution with density function 
given by : 


f(x;a, 3) 


_ r (g + 3 ) g -1 

r(g) r(3) 


(1 - x) 


3-1 


0<x<l , 0<g, 3<“> (26) 


it is most convenient to make use of a relationship between the Beta distribution and 
the F-distribution. That is, if X is a random variable following the F distribution 

with m and n degrees of freedom, then Y = ^ ^ follows the Beta distribution. 

Consequently, ^ 
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Y 


2m 



X. 

1 


2 



2n+2m 



i=2m+l 


(27) 


is the desired Beta random numbers generator, where the X.'s are standard Normal 
random numbers realized by methods in Section III, paragraph A. 

IV. RANDOM NUMBERS FROM DISCRETE DISTRIBUTIONS 


Generation of random numbers from a discrete distribution is handled in a 
manner analogous to that of a continuous distribution. Given a uniform random number 
u and cumulative discrete distribution F(x) , the least value of x for which F(x) >u is 
sought. This X is the desired random number having discrete density function f(x). 


A. Binomial Distribution 
The cumulative binominal distribution is given by 


F(x;n,p) 


i:(^) pNi-p)-'' 

k=0 


(28) 


For a given uniform number u, one may successively evaluate the upper limit x until 
the minimum value of x is found for which 


p’^d - p)"->' iu .2 (29) 

k=0 

An alternate method is based on the sum of n independent Bernoulli random variables. 

By this method, Bernoulli trials, each having probability of success p, are simulated 
with each Bernoulli trial being assigned a one or zero (depending on success or failure). 


2. Because factorials, exponentials, and powers frequently occur in probability 
density functions, care must be taken to avoid computer overflows /underflows 
when computing individual components of the density function. Recursive rela- 
tionships are qmte useful in dealing with these types of problems, but computer 
programs employing recursive relationships are, in general, much slower. 
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Then the sum of these n Bernoulli random numbers will be a random number from a 
Binomial population with parameters n and p. However, it was found that this process 
is relatively slow. 


B . Poisson Distribution 

Poisson random numbers x may be generated from cumulative distribution 
function 


X 

F(x;X) 

k=0 


,k -X 
X e 

k! 


X > 0 


(30) 


by using the same technique which was used for the binominal distribution. An 
alternate method for obtaining Poisson random numbers is by generating uniform 
random numbers u. until the inequality 


k+1 

IT 

i=l 


u. < 
1 


-X 


(31) 


is satisfied, which gives k as the desired Poisson random number. However, this 
method was found to be relatively slow. 


C. Neyman Type-A and Thomas Distributions 

These two distributions belong to the category of cluster (self- exciting) point 
processes and have found application in aerospace engineering, ecology, reliability, 
and forestry. Cluster processes are characterized by a primary (mother) process 
which generates at each point secondary (daughter) events. When only daughter 
events appear in the final process and when primary and secondary distributions are 
both Poisson, the resulting distribution is known as the Neyman type- A counting 
distribution with probability function [8] 


p(n,o.,S) = 


n -am 
e 


s'" 

m! 


m=0 


n = 0,1,2,. .. 
and a, 3 > 0 


(32) 


The parameter 3 represents the rate at which primary events occur and a is the 
average number of secondary counts per primary. The mean and variance is a3 and 
(1 + a)a3, respectively. The Thomas distribution is similar, except that primary 
events are counted in the final process. The probability function is given by 
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n - 0 , 1 , 2 , . . . 


(33) 


n 


p(n;a, 3 ) 


m=0 


(am) 
(n-m)! 


n-m 


m 


3 ill 

e 


m! 


2 

with mean (1 + a ) 3 and variance (1 + 3a + a ) 3 - 

APL programs are provided in Appendix B which generate random numbers from 
these distributions. 


V. CONCLUSION 


Because large scale simulations often require a vast number of random numbers 
from various distributions , emphasis should be placed on speed and accuracy of the 
methods used. We have found the congruential uniform generators presented in 
this report to be both suitable and practical. However, not all possibilities have been 
exhausted and the literature is found to be extensive in this area. Random numbers 
can be generated from nonuniform distributions by finding the inverse to the cumula- 
tive distribution functions in accordance with the probability integral transformation. 
However, when no convenient inverse exists, numerical approximation methods or 
statistical relationships may be used to generate the desired random numbers. Again, 
speed and accuracy should be key factors in choosing between candidate methods. 

Some useful approximations to the eumiolative normal distribution are given in Appendix 
A and APL programs for most of the generators are listed in Appendix B. 
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APPENDIX A 


APPROXIMATION TO THE CUMULATIVE NORMAL DISTRIBUTION 


Because of the extensive application of the Normal distribution in aerospace 
simulation studies, we include here three methods for approximating the cumulative 
standard Normal distribution. Other useful formulae may be found in Reference 9. 
Denote the standard normal density function by 


f(x) = (2^) 


- 1/2 



/2 


and its cumulative distribution function by F(x). Then F(x) may be represented by 
each of the following: 


1 

2 



oo 


E 

n=0 


(_l)nx2n+l 

2%! (2n+l) 


-OO < X < oo 


(A-1) 


1 

1 

2 


-f(x) 

+f(x) 


1_ L_ 2_ 1_ 

x+ x+ x+ x+ 


X 

1 - 


2 o 2 

X 2x 
3+ 5- 




X > 0 


X > 0 


(A- 2) 
(A- 3) 


The last two formulae are called continued fractions, where the conventional notation 


bi+ b2+ •• 


bi+a2 


^2’*’^ 3 
^3^ 


is used to conserve space. Since each formula is exact, we investigate their rates of 
convergence to F(x) for computational convenience. 

Equation A- 2 converges most rapidly for x > 3 while A-1 and A- 3 are preferred 
for 0 < X < 3. 
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APPENDIX B 


APL PROGRAMS 


V R^N CDFI CDF 

[ 1 ] 

[2] L:R^R,+/CDF<UNF 1 

[3] ->-{Q<N^N-l) /L 

[4] ^0 

[5] ft RETURNS N RANDOM NUMBERS FROM THE 

[6] ft DISCRETE DISTRIBUTION FUNCTION 

[7] ft SPECIFIED BY CDF. 

V 

V R^CTRALIM N 

[1] 6+(+/ l + ?ffpl001)^1000;i?««-ff,12 

[2] ^0 

[3] ft THIS FUNCTION GENERATES NORMAL 

[4] ft RANDOM NUMBERS ACCORDING TO 

[5] ft THE CENTRAL LIMIT METHOD. 

1 

V CUMNORM XiCiSiY\AiM\TiZiSO 

[ 1 ] ( 7 -<-( 02 )* 

[2] LO:r^XxJxZx l*L Z-^l+X-Af-<-M+l 

[3] i4'^(2xZ)+ 1+TiA 

[4] S-^S + SO->-SQy--X'x-X'x^{ 1 + 2XW) ^Mx2xl + 2xM 

[5] 7-^A:+z^y 

[6] ■*(M<K)/L0 

[7] P^O 0 0 

[8] P[l]-<-l-(7x(*-yxyf2)^y 

[9] P[2]-«-0.5+C7x(*-yxy^2)xy^>i 

[10] P[3]^0.5+Cx5 

[ 11 ] -»-0 

[12] ft THREE FORMULAE FROM THE HANDBOOK OF MATHEMATICAL 

[13] ft FUNCTIONS, EDITED BY ABROMOWITZ, 1970 , ARE PROVIDED. 

[14] ft THE FIRST TWO ARE CONTINUED FRACTIONS AND THE THIRD IS 

[15] ft A POWER SERIES. P[l] {EQ. 26.2.14, PAGE 932 , GIVES 

[16] ft P(y), y>0. P[2] ALSO GIVES P(X) AND IS EQ. 26.2.15. 

[17] ft P[3] IS EQ. 26.2.10. [1] CONVERGES MOST RAPIDLY FOR 

[18] ft y>3 WHILE P[2] AND P[3] ARE PREFERRED FOR 0<X<3. 

[19] ft PROGRAMMED BY L. HOWELL, 19 MAY 82. 

7 

V E^PARAM GENLOGNORM N ’,MU %V AR 

[1] VAR^PARAMl2'\',MU->^PARAMil'S 

[2] E^*MU+{VAR*0 .S)^NORM N 

[ 3 ] ->-0 

[4] ft GENERATES N LOG-NORMAL RANDOM 

[5] ft NUMBERS WITH PARAMETERS MU AND 

[6] ft VAR {MU AND VAR ARE THE MEAN 

[7] ft AND VARIANCE OF THE DEFINING NORMAL 

[8] ft DISTRIBUTION) . 

7 
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V E->rEASTING N iC ; D iUO iNU i DE ;T 

Cl] C-<-2. 515517 0.802853 0 . 010328 ;P'^1 1.432788 0.189269 0.001308 

[2] I/0-^"0. 5001+(.9/|^plO00)*1000 

[3] r-^(-2x«|i/0)*0.5 

[4] NV->-G[:il+T^Cl22+T^CldJ 

[5] Z)ff-<-Z?[l]+2’xZ?C2]+2’xZ)[3]+rx2?[4] 

[6] R^(>^U0)^(T-IIU*DE) 

[7] -►O 

[8] fi GENERATES NORMAL RANDOM NUMBERS 

[9] fi USING HASTING^ S RATIONAL APPROXIMATIONS 

V 


V CD^Ftlf NJilYMAJS A^ZitfjMiXiXsN;NPiRiS 
[!•] NP^*-V^l-*-A 

[2] N->-0 

[3] Ll:N^N+l;S'>-i0iM->-0 

[4] L2:X-<-liM*-l+M 
C5] II->-l 

[ 6 ] L3:X-^X^AyMiN+l-II 

[7] ^(N^II-*-II+l) /L3 

[8] X-^X^*-A^M 

[9] y^-iiJ^l 

CIO] L4;y-fJxAr*Af+i-j 
Cll] -^iM^J^J+l) /L^ 

Cl2] Y^y^*-w 
Cl3] 5-*-5.yxy 

Cl4] -»•( (i?<0)v(i?>lff~12)v(A/S5) )/I2;/?>^-/"2+5 

Cl5] NP-*-NP^ + /S 

Cl6] -*-(0.999> + /i^P)/Ll 

Cl7] CDP^+\NP 

Cl8] -►O 

Cl9] n COMPUTES THE CUMULATIVE NEIMAN-TIPE-A 

C20] 4 PROBABILITY DISTRIBUTION FUNCTION. 

C21] Pi W IS AVERAGE NO. OF PRIMARY (.MOTHER) EVENTS 

C22] Pi A IS AVERAGE NO. OF SECONDARY {DAUGHTER) EVENTS PER PRIMARY 

C23] n PROGRAMMED BY LEONARD HOWELL » AUGUST 4. 1982. 

7 


V R^NORM S 

Cl] P-^"6 + 0. 0001 x + /?5p 10000 ^2. 

7 
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V PDF^PDFFRCDF CDFiY 

[ 1 ] PDF->-CDF’-YiY’^0,l^~l<^CDF 

V 

V CDF-^POISSON STiPQiPiK 

Cl] CDF-*-PQ-<-*-ST 

C2] K'^P-^l 

[ 3 ] LOzP^P^STiK 

C4] CDF->-CDFA~l^CDF)+PxPO 

C5] K'>-K+l 

[6] ^(0.9999>'l + CZ?F)/Ii0 

[ 7 ] -►0 

[ 8 ] ft RETURNS THE CUMULATIVE POISSON 

[9] ft DISTRIBUTION FUNCTION WITH 

[10] ft PARAMETER ST. 

7 

7 R^MU POISl KiFiX 

[1] R^{*~MU)^(MU*X)ilX^0,\K 

C2] F->-+\R 

[3] ^0 

[4] ft CUM. VIST. FUNCT. FOR POISSON 
7 

7 R^MU P0IS2 KiFUPiX 

[1] X^0;P-^R->-i*-MU) i 

[2] L:P^MU^P*X^X+1 

[ 3 ] R^R,P 

[4] ■^(K*X)/L 

[5] F1^+\R 
7 


7 R^POLAR N 

[1] R-^(?R)iR->-(2j x/i|?,0.5)pl0737418 24 

[2] R*-Npil 2o.O(02)x"l 0 + i?)x("2 "2».x®l 0+i?)*0.5 

[ 3 ] ■►0 

[4] ft GENERATES RANDOM NUMBERS USING 

[5] ft THE BOX-MUELLER METHOD, OFTEN 

[ 6 ] ft REFERRED TO AS THE POLAR METHOD. 

7 

7 R^NP RBINOl NiAiI;UiK 

[ 1 ] 

[2] ii0:i4+-0;J-^0 

[3] L:i/-<-{"l + ?1001)*1000 

[4] K^U^NPl22 

[5] A^A+K 

[6] -*-(il?P[l]iI-^I+l)/i 

[7] R-*-R,A 

[ 8 ] -*-(.0*N-^N-l)/L0 

[9] -»-0 
7 
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V E^flP RBIN02 NiFiMiP;PFiU iXiY 

[1] X<-0 ,iMiM^NPli:ilP^NPL2l 

[2] PF^iXlM)yi{P*X)>^(l-P)*M-X 

[3] F*-+\PF 

[4] i/-«-("l + ?ffpl001)vl000 

[5] R-^+/U°.>F 

[6] ->0 

[7] n GENERATES BINOMIAL RANDOM NUMBERS 

[8] fi USING THE INVERSE METHOD. 

[9] n ffP[l] IS THE SAMPLE SIZE AND 

[10] fl NPL2) IS THE PROBABILITY 

V 

V R-<-NP RBIN03 NiM;PiU 

[ 1 ] M<-NPLi:\iP^NPL2l 

[2] i/-«-(~l + ?(ff,Af)pl001)TlOOO 

[3] R^+/U>P 

[4] ->0 

[5] n GENERATES BINOMIAL RANDOM NUMBERS 

[6] f\ AS THE SUM OF BERNOULLI TRIALS. 

[7] R ffPCl] IS THE SAMPLE SIZE AND 

[8] fl ffp[2] IS THE PROBABILITY . 

V 


V R^MU RP0IS3 NiAiliKiU 

[1] 

[2] LOzA^liK^O 

[3] LzU^i 1+?1001 ) vlOOO 

[4] A*-A^U 

[5] K^K+1 

[6] -^iA>(*-MU)) / L 

[7] R-^RA 1 + K) 

[8] -^iN>I-<-I + l) / LO 

[9] ^-0 

[10] fl GENERATES POISSON RANDOM NUMBERS 

[11] fl BY TAKING THE PRODUCT OF UNFORM 

[12] fl RANDOM NUMBERS AND CHECKING FOR 

[13] fl FOR THE INEQUALITY EXP(,-MEAN) 

V 

V R^MU PP0JS4 N iX’,AiP-,FiU iM 

[ 1 ] X-<rOiA^P*-{*-MU) 

[2] L:A^MU^AiX^X+l 

[ 3 ] P-<-P , A 

[4] -^(XxfM^MU+HyiMUi^O .5) / L 

[5] P-«-P,>l 

[6] F^+\P 

[7] P-<-("l + ?ffpl001) vlOOO 

[8] R^+/Uo.>F 

[ 9 ] ->-0 

[10] fl USES THE INVERSE CUMULATIVE 

[11] fl METHOD FOR DISCRETE RANDOM 

[12] fl NUMBERS 



ISTDUNFSUMWl 

7 STDUNFSUM It ;K ;X ; S i Z iPRT lY lUX i SIG ;X1 lU 

[1] UX->-Nr2iSIG^iNil2)*0 .5 

[2] X^-UXiSIG;K^~l+i{N+l) iY^\OiPRT->-iO 

[3] L:Z^iXl-K)^U^Xl>KiXl^UX+SIGxX 

[4] S*-SIG^(.-/ (KlN)y<Z*N-l)^liN-l) 

C5] Y^Y,S 

[ 6 ] PRT^PRT , X 

[7] A:^J+0.0 5 

[8] ^(XiUXiSIG) /L 

[9] 0 DRAW Y VS PRT 

[10] ■►o 

[11] fi COMPUTES THE PROBABILITY DENSITY 

[12] fi FUNCTION OF THE STANDARDIZED 

[13] fi SUM OF N INDEPENDENT UNIFORM 

[14] n RANDOM VARIABLES BY TAKING THE 

[15] ft INVERSE OF ITS MOMENT GEN- 

[16] ft ERATING FUNCTION. PROGRAMMED 

[17] ft BY LEONARD HOWELL 08/04/82. 

7 


VSTUDENTWV 

7 T^V STUDENT P \C \D \R\XP \G\ \G2 \GZ \Gi\ \G iVP \DUD2 

[1] C*^2. 515517 0.802853 0 . 0 1 0 3 2 8 ; 1.432788 0.189269 0.001308 

[2] i?*^("2x®P)*0 . 5 ;7P-^7*“l ~2 3 "4 

[3] XP->-R-Dl^D2iDl-<r+/C^R*0 1 2 ; P2-<-+/Pxp*o 12 3 

[4] G1*h( + /ZP*3 l)f4 

[5] G2-<-( + /(ZP*5 3 px5 16 3)f96 

[6] G3^(+/3 19 17 15xZP*7 5_3 1)^384 

[7] G4-^( + /79 776 1482 1920 945xZP*9 7 5 3 1 )^92160 

[8] (?-«-Gl ,G2 ,G3 ,G4 

[9] T^XP+ + /G>^VP 

[10] -*-0 

[11] ft THIS FUNCTION APPROXIMATES THE INVERSE STUDENT -T 

[12] ft DISTRIBUTION AND IS USEFUL IN AUTOMATIC LOOK-UP OF 

[13] ft CRITICAL VALUES IN HYPOTHESIS TESTING. FORMULA 

[14] ft 26.7.5, PAGE 949 OF ABRAMOWITZ' S HANDBOQK-Ql 

[15] ft MTQ.MA11Q.AL_FUNQ.ZIQ.IS. TS USED IN CONJUNCTION WITH 

[16] ft FORMULA 26.2.23 {APPROXIMATES THE INVERSE NORMAL 

[17] ft DISTRIBUTION) , PAGE 933. PROGRAMMED BY L. HOWELL, 

[18] ft JUNE 23, 1982. 

7 
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VTHOMASLUl'7 

V CDF^THOMAS I NP ;M -,0 ;E i ED ‘,1 ‘,P iT ;IJ iJl 

[1] M^INPlllPi AV NO OF MOTHER EVENTS {PER MS) 

[2] D^INPl2'\(\ AV NO OF DAUGHTER EVENTS {COUNTS) 

[3] E^{*-M)>^1,M 

[4] r-f-SLl] 

[5] CDF-<rEll'\ ,E\.ll + E\.2'\^ED*-*~D 

[6] 1^2 

[7] 

[ 8 ] E-^EA l^E)xMiI 

[9] L2:P^P^DiJ 

[10] ^{I>J->-J+l) /L2 

[ 11 ] P-*-P^ED 

[12] T->-PxEL2l 

[13] J<-1 

[14] L3:P->-Px{ {JliJ)*IJ)xIJxEDiDxJl-<-J + l;IJ->-I-J 

[15] T-<rT + PxELJ+2l 

[16] ■^{I>J-^J+1)/L3 

[17] CDF^CDFA 1\CDF)^T 

[18] I^J+1 

[19] ->-( 0 . 9999>"1'K;2?P) /PI 

[20] ^-0 

[21] PI RETURNS THE CUMULATIVE THOMAS DISTRIBUTION 

V 


SlUNFimN 

V U^UNF S;N 

[1] U^iO;N<-x/S 

[ 2 ] L: U^U , UNFDxUNFX^UNFM \ UNFRxUNFX 

[3] -y{0<N^N-l)J L 

[4] U^SqU 

[ 5 ] ->-0 

[6] PI STANDARD UNIFORM NUMBER GENERATOR 

V 


NUNFICLUlN 
V UNFIC SEED 

[1] Z7i7PM-^214748 3 647 

[2] f/i7Pi?^29 9 0 39 47 

[3] UNFX^SEED 

[4] UNFD^iUNFM 

[5] ^'0 

[6] PI INITIALIZES PARAMETERS FOR THE 

[7] PI STANDARD UNIFORM GENERATOR. 


V 



VUNFSUMIQIV 

V UNFSUM_NiKiXiSiZiUiPRTiY 

[1] X-<-OiK-<- 1+\{N+1) iY-<riO;PET*-iO 

[2] L:Z-<^(X-K)xU-<rX>K 

[3] (KlN)xZ*N~l)il (N-1) 

[4] 

[5] PRT‘>rPRT,X 

[ 6 ] 

[7] -^{XiN)/L 

[8] 0 DRAW Y VS PRT-N^2 

[9] fl COMPUTES INVERSE OF THE MOMENT 

[10] fl GENERATING FUNCTION OF THE SUM OF 

[11] PI INDEPENDENT UNIFORM VARIABLES 

V 


VUNIFORMlW]! 

V UNIFORM N\I%U',A-,S 

[1] R^\0-,I^l-,U-<-SEEDiA^l%Q01 iMOD*- 1 + 2*31 

[2] L:U^M0D\AxU 

[3] SHMxU^MOD 
[ 4 ] R^R , S 

[5] -^(N>I-^I+1) /L 

[ 6 ] +-0 

[7] ft STANDARD UNIFORM NUMBER GENERATOR. 


V 



REFERENCES 


1. Lehmer, D. H.: Proc. 2nd Symposium on Large-Scale Digital Computing 

Machinery. Cambridge, Harvard University Press, 1948, pp. 141-146. 

2. Marsaglia, G. : The Structure of Linear Congruential Sequences; Applications 

of Number Theory to Numerical Analysis. S. K. Zaremba (ed.), Academic 
Press, 1972. 

3. Atkinson, A. C. : Test of Pseudo-random Numbers. Applied Statistics, Vol. 

29, No. 2, pp. 164-171. 

4. Carroll, Stan N.: Aerospace Engineer, Dynamics and Trajectory Analysis 

Branch, Control Systems Div. , Systems D 3 mamics Lab., MSFC , NASA, Huntsville, 

AL, personal communication, 1982. 

5. Harris, B.: Theory of Probability, Addison-Wesley , 1966. 

6. Hastings, C.: Approximations for Digital Computers. Princeton University 

Press, Princeton, N.J., 1955. 

7. Box, G. E. P. and Mueller, M. E.: A Note on the Generation of Random Normal 

Deviates. Ann. Math. Statist., Vol. 29, pp. 610-611. 

8. Neyman, J. and Scott, E. L. : Processes of Clustering and Applications, in 

Stochastic Point Processes: Statistical Analysis, Theory, and Applications, P.A. 

W. Lewis (ed.), Wiley- Interscience , New York, 1972, pp. 646-681. 

9. Abramowitz, M. and Stegun, 1. A.: Handbook of Mathematical Functions with 

Formulas, Graphs, and Mathematical Tables. National Bureau of Standards, 
Applied Mathematics Series, 55, U. S. Government Printing Office, Washington, 
1964. 


21 



BIBLIOGRAPHY 


Aitchison, J. and Brown, J. A,: The Lognormal Distribution. Cambridge University 

Press, 1963, 

Feller, William: An Introduction to Probability Theory and Its Applications. Vol. 1, 

Wiley, 1964. 

Hahn, G. J. and Shapiro, S. S.: Statistical Models in Engineering. Wiley, 1967. 


22 



1. REPORT NO. 

NASA TP-2105 


2. GOVERNMENT ACCESSION NO. 


4. TITLE AND SUBTITLE 

Generation of Pseudo -Random Numbers 


7- AUTHOR(S) 

Leonard W, Howell and Mario H- Rheiirfurth 


9. PERFORMING ORGANIZATION NAME AND ADDRESS 


George C. Marshall Space Flight Center 
Marshall Space Flight Center, Alabama 35812 


3. RECIPIENT’S CATALOG NO. 

5. REPORT DATE 

December 1982 

6. PERFORMING ORGANIZATION CODE 
8. PERFORMING ORGANIZATION REPORT # 
10. WORK UNIT NO. 

M-595 . 

|1 1. CONTRACT OR GRANT NO. 


12. SPONSORING AGENCY NAME AND ADDRESS 

National Aeronautics andCSpace Administration 
Washington, D.C. 20546 


|l3. TYPE OF REPORT & PERIOD COVERED 

Technical Paper 


14. SPONSORING AGENCY CODE 


15. SUPPLEMENTARY NOTES 

Prepared by Systems Dynamics Laboratory, Science and Engineering Directorate. 


16. ABSTRACT 

The generation of pseudo-random numbers from specified probability distribu- 
tions has found extensive application in Monte Carlo simulations. Because these 
simulations often require a large number of calculations, the time required to 
generate pseudo-random numbers has become a major factor. At the same time it 
is essential to guarantee accuracy and statistical randomness of the sequence of 
generated numbers. This report provides practical methods for generating 
acceptable random numbers from a variety of probability distributions which are 
frequently encoxHitered in engineering applications. 


17. KEY WORDS 

Random Number Generators 
Monte Carlo Simulations 
Pseudo-Random Numbers 
Congruential Sequences 
Uniform Random Numbers 


19. SECURITY CLASSIF. (of this r*port) 

Unclassified 


18. DISTRIBUTION STATEMENT 

Unclassified - Unlimited 
STAR Cagetory: 65 


20. SECURITY CLASSIF. (of this psgs) 


21. NO. OF PAGES 

22. PRICE 

Unclassified 


27 

A03 


For sale by National Technical Information Service, Springfield. Virginia 221 Gl 


NASA-Langley, 1982 



National Aeronautics and 
Space Administration 


THIRD-CLASS BULK RATE 


Washington, D.C. 
20546 

Official Business 

Penalty for Private Use, $300 





POSTMASTER: If Und el iver able (Section 158 

Postal Manual) Do Not Return 




